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ABSTRACT 

The measurement error of pulse times of arrival (TOAs) in the high S/N limit is 
dominated by the quasi-random variation of a pulsar’s emission profile from rotation 
to rotation. Like measurement noise, this noise is only reduced as the square root 
of observing time, posing a major challenge to future pulsar timing campaigns with 
large aperture telescopes, e.g. the Five-hundred-metre Aperture Spherical Telescope 
and the Square Kilometre Array. 

We propose a new method of pulsar timing that attempts to approximate the 
pulse-to-pulse variability with a small family of ‘basis’ pulses. If pulsar data are in¬ 
tegrated over many rotations, this basis can be used to measure sub-pulse structure. 
Or, if high-time resolution data are available, the basis can be used to ‘tag’ single 
pulses and produce an optimal timing template. With realistic simulations, we show 
that these applications can dramatically reduce the effect of pulse-to-pulse variability 
on TOAs. Using high-time resolution data taken from the bright PSR J0835—4510 
(Vela), we demonstrate a 25-40% improvement in TOA precision. Crucially for pulsar 
timing applications, we further establish that these techniques produce TOAs with 
gaussian residuals. 

Improvements of this level halve the telescope time required to reach a desired 
TOA precision. Although some gains can be achieved with existing data, the greatest 
improvements result from the ‘tagging’ approach, which in turn requires online or 
posthoc analysis of single pulses, an important consideration for the design of future 
inst rument at ion. 
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1 INTRODUCTION 


From the discovery of the very first pulsar by |Hewish et al.| 
119681, two key features of pulsar emission were evident: (1) 


the extraordinary regularity of the pulse train and (2) the 
wide variety of the shapes and intensities of the single pulses 
of the train. The former admits a number of important ex¬ 
periments via pulsar timing. E.g., the timing of unrecycled 
and mildly recycled pulsars provided the first evidence of 
gravitational radiation (Taylor, Fowler fc McCulloch||1979 1 


and other stringent tests of general relativity ( Kramer fc 
WexpOdg Ransom et al.|[2014 l, as well as the discovery of 


magnetospheric reconfiguration (Kramer et al. 2006 Lyne 


et al.|20l0 Hermsen et al.|2013 1. Timing millisecond pulsars 


(MSPs) probes physics at supra-nuclear density via neutron 
star mass measurements | Lattimer fc PrakashpOOl Demon 


est et al.|2010 Antoniadis et al.|2013 l and enables searches 

for low frequency gravitational waves (e.g. Yardley et al. 


2011 Demorest et al.|2013 Shannon et al.|2013l. 

The precision of timing experiments is limited by a va- 
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riety of stochastic noise phenomena. The effects of propa¬ 
gation through the interstellar medium are strong but are 
largely correctable with broadband observations; see e.g. 
Keith et al. (20131. For truly faint pulsars, the limiting 


factor is the white radiometer noise added to the observed 
pulse profile. For most pulsars, however, the primary limit 
arises from ‘timing noise’, a poorly-understood collection of 
phenomena in which the residuals of pulse times of arrival 
(TOAs) show time-correlated structure, i.e. red noise. Be¬ 
cause timing noise seems to be correlated with spin-down lu¬ 


minosity and other proxies for pulsar age (e.g. Hobbs, Lyne 


& Kramer 20101, it most strongly affects unrecycled pul¬ 


sars. However, its presence has also been detected in MSPs 
( Kaspi, Taylor fc Ryba|1994 Manchester et al.|2013 ), where 
it dilutes sensitivity to the correlated residuals expected to 
be induced by gravitational waves. 

A second phenomenon—which in the literature has 
been called both ‘jitter noise’ ( [Shannon fc Cordes||20i0| 
and ‘stochastic wide-band impulse modulated self-noise’ 
(SWIMS; jOslowski et~aL][20TI| — arises from the second 
fundamental property of pulsar emission, its pulse-to-pulse 
variability. Despite the formidable nomenclature, the phe- 
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nomenon is simple: because we only record a finite number 
of pulses in an observation, the resulting profile is randomly 
distorted from the mean, biasing the measurement of phase. 

Like radiometer noise, which fundamentally limits any 
radio experiment, jitter/SWIMS seems to primarily induce 
white noise in timing r esiduals and, for a given observation , 
decreases as oc 1/Vi (Shannon, Oslowski & et al. 20141. 
Thus, the ratio of jitter/SWIMS to radiometer noise, (jj/ar, 
depends only on the properties of the pulsar and the receiv¬ 
ing system, and historically, the jitter/SWIMS phenomemon 
has only been important for the brightest pulsars. But in¬ 
creasingly sensitive receiver systems have brought the issue 
to the forefront and several studies have been published in 
recent years. In particular, Oslowski et al. I® present 
both an excellent overview of the phenomenon (Dear Inter¬ 
ested Reader: run, do not walk, to print out a copy of this ex¬ 
cellent paper) and a study of its effect on PSR J0437—4715, 
by two orders of magnitude the brightest MSP. Those au¬ 
thors conclude jitter/SWIMS dominates measurement un¬ 
certainty, contributing at about thrice the level of radiome¬ 
ter noise. 

Because PSR J0437—4715 is so bright, it is a bellwether 
for future timing experiments with large aperture telescopes, 
e.g. FAST (the Five-hundred-metre Aperture Spherical Tele¬ 
scope) and the SKA (Square Kilometre Array). Many pul¬ 
sars whose timing is currently limited by telescope sensi¬ 
tivity will be limited by jitter/SWIMS, with severe impli¬ 


cations for high precision timing experiments (e.g. Hobbs 


et al. 20141. In particular, without the ability to mitigate 


jitter/SWIMS, the only tractable solution will be to time 
large numbers of MSPs at lower precision. This, in turn, 
places substantial pressure on the design of FAST and SKA, 
viz. slew time for the former and the availability of many tied 
subarrays for the latter. 

In this work, we consider a new method for treating jit¬ 
ter/SWIMS based on analysis of single pulses. By approxi¬ 
mating the pulse-to-pulse variability as arising from a finite 
family of pulse shapes, we attempt to estimate the particu¬ 
lar realization of sub-pulse structure in a given observation 
and from this produce a more accurate TOA. We give an 
overview of this approach in the context of pulsar timing 
in In the following we briefly describe high-time res¬ 
olution observations of the Vela pulsar (PSR J0835—4510) 
before going on to discuss the approximation of its pulse-to- 
pulse variabilty by construction of a template basis in 
In j|5.1| and j]5.2| we discuss application of the basis method 
to simulated and real data collected both with and with¬ 
out single pulses, respectively, and demonstrate substantial 
improvement in TOA precision in both cases. Finally, we 
summarize our results and discuss their implication for and 
possible implementation in future timing experimets. 


2 PULSAR TIMING 

Elementary pulsar timing comprises a series of measure¬ 
ments of the time of arrival of a pulse in the laboratory 
frame. Rather than working directly with time series, it is 
typical to work instead in terms of pulsar phase and to mea¬ 
sure the phase offset, 5, between an observed pulse shape, p, 
and a template pulse shape, t. If pulse shape and template 
are realized as A/,-bin vectors p and t, then the probability 


density fnnction (pdf) for p is given by 


/(pi 5,s,t,(j)= 


1 = 1 ^ 


exp 


Pi - sti{d) 




( 1 ) 

where a is the white radiometer noise and both data and 
template are baseline-subtracted. Implicit in this formula¬ 
tion is the assumption that pulsar emission can be well- 
described as amplitude-modulated noise (Rickett 19751. 
Most pulsars raise the system temperature negligibly, and 
a can be taken as a constant independent of phase and es¬ 
timated from the baseline variance or receiver parameters, 
and this is done in ‘classical’ pulsar timing. For bright pul¬ 
sars, especially those we consider here, the peak intensity 
may exceed the system equivalent flux density, and a is in¬ 
creased by pulsar ‘self-noise’ (e.g. Kulkarnl||1989 Gwinn & 


Johnson|2011 1, becoming dependent on phase and unknown 


a priori. As usual, once an observation provides a fixed real¬ 
ization of p, / can be viewed as the likelihood function for 
the parameter(s), T((5), and an estimator, 5 determined via 
maximum likelihood. In pulsar timing, the likelihood func¬ 


tion is typically evaluated in the Fourier domain (Taylor 


19921. 


In the above formulation, the pdf / is conditioned on 
the template t, so if our assumptions about t are wrong, 
S will be biased. In pulsar timing, t is typically taken to 
be the mean pulse profile obtained from a long integra¬ 
tion, or a smoothed/analytic approximation thereof. Be¬ 
cause t changes with each pulse, the observed profile p only 
approaches the mean asymptotically, with fractional error 
oc 1/Vi- The increased scatter—^jitter/SWIMS—is a mani¬ 
festation of the bias incurred by using a time-averaged tem¬ 
plate. 

There are two fundamentally different paths to dealing 
with this uncertainty, both involving a modification of the 
likelihood. As presented by Oslowski et al. (20111, one ap¬ 
proach is to treat the pulse-to-pulse variability as correlated 
noise. This solution is elegant as it captures both the in¬ 
creased variance from phase-dependent system temperature 
(their ‘Regime 2’) and the pulse-to-pnlse variability (their 
‘Regime’ 3). It comes at the cost of a more complicated like¬ 
lihood, as the phase bins are no longer independent: 


f{p\S,s,t,a) 


1 

, exp 

VStt det C 



•Cr 


( 2 ) 


where C is the non-diagonal covariance matrix of the pro¬ 
file bins and r are the residuals to the template. Gener¬ 
ally, C must be measured from the data. [Oslowski et al.| 
( |2011[ ) measured C and used a principle component analysis 
(PGA) to recover some of the jitter/SWIMS distortion and 
improve the timing precision of PSR J0437—4715 by about 
20%. [Oslowski et al.| ( |20f^ extended the method to include 
polarimetry and achieved a 40% improvement in precision. 

The second approach, which we adopt, views the tem¬ 
plate itself as an unknown to be determined from the data. 
The latter remain uncorrelated, though the errors may be 
heteroscedastic due to self-noise. Formally, t is a nuisance 
parameter in Hilbert space, and even approximating it as t 
incurs many additional degrees of freedom. Moreover, al¬ 
though the pulse-to-pulse variability of some pulsars has 
been stndied in detail, it remains a poorly understood, 
stochastic process. 
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In this work, we make the simplifying assumption that 
a given single pulse is drawn from a hnite family of Nb pulse 
shapes, {6}. For an integration over many pulses, this is 
equivalent to expanding t in a small number of basis func¬ 
tions: 

JVt 

(3) 

where the Si will be determined by the number and fluxes 
of each pulse family. 

The focus of the work below is to test this approach in 
two scenarios. In both cases, we have a ‘training’ set of single 
pulse data from which a template basis can be constructed. 
(We discuss one approach to basis construction in SQ) In the 
hrst scenario, all other timing observations are carried out in 
‘fold mode’, in which many single pulses are coadded. This 
configuration is typical of pulsar timing observations. In this 
case, we attempt to reconstruct the single pulse information 
by estimating the basis components of equation from the 
fold-mode data. 

In the second scenario, we assume we have at least 
limited access to single pulse profiles for each timing ob¬ 
servation. (We discuss precisely what information might be 
needed in (j^) In this case, we can build up an optimal pro- 
Hle for each observation by mapping the single pulses to our 
basis functions, and we find that such an optimal profile 
improves the timing precision even further. 

To test this method, we use high time resolution obser¬ 
vations of Vela, PSR J0835—4510. Vela is the brightest pul¬ 
sar in the sky, and like PSR J0437—4715 is an excellent test 
case for probing the jitter/SWIMS limit on current receiver 
systems. We describe the observations and data reduction 
i n jj3 Both cases require a template basis, and in §4.1| and 
$ |5.3 we discuss two effective methods for its creation. 


3 OBSERVATIONS 


On 4 July 2012, we used the 64-m Parkes telescope and the 
CASPSR backend to record baseband voltage samples from 
the two orthogonal, linearly-polarized feeds of a receiver 
tuned to 2820 MHz. The 800 MHz sampling rate provided 
useful bandwidth of about 320 MHz after filtering. Data were 
recorded in 8-second blocks, and although the full observa¬ 
tion lasted about two hours, many blocks of data were lost 
from the buffer due to the failure of a disk in the on-line 
recording system. Ultimately, we obtained about 30 min¬ 
utes of data, or 21,184 pulses, in non-contiguous 8-second 
blocks (see Figure]^. We note that 8 seconds is substan¬ 
tially longer than any reported pulse-to-pulse correlations 
(e.g. Krishnamohan &: Downs|1983 |. 

Using DSPSR ( [van Straten k, Bailes||20ll| , we synthe¬ 
sized coherently dededispersed filterbanks of the Stokes pa¬ 
rameters of each single pulse with spectral and time res¬ 
olution of 512 channels and 1024 phase bins (~87 qs), re¬ 
spectively. We corrected the filterbanks for the complex dif¬ 
ferential gain between the two feeds and for the bandpass 
response using observations of a pulsed noise source, and we 
hnally integrated the filterbanks over frequency and polar¬ 
ization to produce monochromatic total intensity (stokes I) 
profiles. The S/N of the resulting single-pulse profiles appear 
in Fig. Vela is truly outrageously bright. 
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Figure 1. The S/N of our single-pulse observations. The peak 
S/N for a single phase bin appears to the left in salmon, while 
the total S/N summed over the onpulse region appears to the 
right in black. The bottom panel shows the time elapsed during 
the observation and missing intervals; the top panel shows the 
cumulative S/N distribution. 


The precise topocentric pulsar frequency at the time of 
observation is unknown a priori. We measured it by extract¬ 
ing a TOA from each pulse and htting a linear model. We 
computed the phase offset of each TOA to this linear model 
and then rotated the pulse profiles by this offset, equivalent 
to folding the data at the exact pulsar frequency. It is im¬ 
portant that this phase alignment matches the true pulsar 
rotational phase with adequate precision. Because we have 
used more than 10"* pulses to estimate it, the phase align¬ 
ment of any individual TOA will be better by a factor of 
100 than the TOA measurement uncertainty. As we shall 
see below, this is adequate for our analysis. 

We monitored the gain stability during the observation 
by analyzing the baseline (offpulse) noise in our single pulse 
profiles, dividing them into 276 (748) onpulse (offpulse) bins. 
The mean offpulse level drifted by about 5% over the obser¬ 
vation, indicating the absolute gain was stable to at least this 
level. Indeed, it was likely more stable, as we found that the 
ratio of the mean to the standard deviation, which is linearly 
proportional to the system gain, k = 4.2282(7) x 10“^, did 
not vary over the observation. We thus expect S/N to be a 
good proxy for flux density. We note our observed k is also in 
agreement band-limited noise, viz. the radiometer relation, 
which predicts k = 4.23 x 10“^ for 320 MHz of bandwidth. 

The dual-band 10/50 cm ( Granet et al.|2005 l observing 
system we used has a system equivalent flux density (SEFD) 
of 25-30 Jy over the observing band, and the peak intensity 
of Vela reaches about 9 Jy (mean about 220 mjy). Thus, for 
a typical pulse, the system noise level is only increased by 
about 5% at peak, a negligible effect. However, peak inten¬ 


sity of the brightest pulses (giant micropulses, see Johnston 


et al. 20011 can be substantially brighter and exceed the 


SEFD, increasing the noise by more than 40%. We address 
this increased noise level further in 

Two important systematic effects could affect a single 
pulse analysis. First, broadband RFI can mimic sub-pulse 
structure. Such RFI is easily identified for PSR J0835—4510 
by its lack of dispersion, and moreover is largely absent from 
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Figure 2. An example template basis. The leftmost basis function is constructed from pulses most resembling the mean, while the 
remainder are based on exceptional single pulses. The salmon, vertical dashed lines mark the same fiducial phase, viz. the maximum 
intensity in the mean profile. Several of the basis functions represent pulses with very bright sub-pulses at particular phases, while others 
capture unusual mixtures of components. 


the 10 cm band; we identified and discarded only four single 
pulses affected by impulsive RFI. Second, as the antenna 
tracks the source, the parallactic angle relative to the feed 
symmetry axis evolves, changing the illumination of the two 
feeds. Correcting such effects requires sophisticated model¬ 
ing of the receiver response (e.g. |van Straten||200'^ |2013[ ). 
But, as the receiver is co-axial with little cross-polarization 
and the range of parallactic angles covered during the ob¬ 
servation is small, we expect this effect to be minimal. 


4 A ‘PULSE SHAPE’ TEMPLATE BASIS 


Here, we discuss an approach for constructing a template ba¬ 
sis based primarily on capturing the various observed pulse 
shapes. In |5.3| we discuss an approach based on the flux 
density of single pulses. In general, we expect any such a ba¬ 
sis to only encapsulate a fraction of the true pulse-to-pulse 
variability. If single pulses are actually formed from random 


realizations of sub-pulses or ‘shots’ (Cordes 1975 Hankins 
|et al.|[2003| ), there is no expectation that a family of sin¬ 
gle pulse shapes should emerge. Moreover, rare but bright 
pulses, such as the giant micropulses of PSR J0835—4510 
or the true giant pulses of PSR B1937-I-21 ([Cognard et al. 


1996), substantially affect TOA measurements but cannot 


be faithfully represented by a small family of basis func¬ 
tions. 

Nonetheless, the question is not one of perfection, but 
of efficacy, and below we show that a modest basis set can 
capture a substantial portion of single pulse variability. 


4.1 Construction 

An effective template basis is one that reproduces the pulse- 
to-pulse variation reasonably well without admitting too 
many degrees of freedom. In principle, for a fixed number 
of basis functions and a given set of single pulses, one could 
determine an optimal basis, i.e. one that minimizes the mean 
squared error between each single-pulse and the best-fitting 


basis function. However, this is an extremely challenging op¬ 
timization problem! 

Instead, we develop a basis iteratively. The initial basis 
comprises the mean template, and we begin by identifying 
the pulse least similar to the mean (see below). We then 
divide the single pulses into those more similar to the mean 
and those more similar to the outlier, and average those 
pulses to produce a basis with two templates. Next, we select 
the pulse least like either of these templates, re-classify the 
single pulses, produce a basis with three templates, and so 
on and so forth. 

Generally, single pulse flux densities follow a lognormal 
distribution ( Burke-Spolaor et al.|2012 I, while some pulsars 
displaying giant pulses and micropulses follow power law 


statistics, (see, e.g., Cairns, Johnston & Das 2001 Cairns 


2004). To avoid constructing our basis from rare bright 


pulses, which would naturally be the strongest outliers due 
to their high S/N, we adopt a ‘reduced’ x^) 


Np Np 

= (4) 

4 = 1 4=1 


which measures the distortion between the observed pulse 
shape and the template normalized to the S/N. 

We note that this procedure does not necessarily con¬ 
verge to a unique basis. With each iteration, the basis func¬ 
tions are revised, and profiles ‘belonging’ to one basis func¬ 
tion may be shuffled to others during the next classification. 
In practice, we find the the functions and classifications set¬ 
tle to a stable result after a few iterations. 

To construct a mean pulse shape template, we co-added 
the first 10"^ Vela pulses. We then applied the basis con¬ 
struction algorithm outlined here to the same pulses to com¬ 
pute template bases with varying numbers of basis functions. 
Templates from an example resulting basis appear in Figure 
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Figure 3. The total mean squared error (x^)i normalized by the 
total degrees of freedom, for the 10^ single pulses as a function 
of the number of functions in the template basis. The black solid 
line traces shows the ‘pulse shape’ basis (ij^ and the scales 
approximately as The salmon, dashed line shows the same 

trend for the ‘flux density’ basis of j]5.3| 


4.2 Application and Efficacy 


Here, we explore how well the basis encapsulates the pulse- 
to-pulse variability. First, we analyze the reduction in total 
mean squared error (x^) of the set of single pulses relative to 
the mean template. For each pulse, we compute relative 
to either the mean template or the best-fitting basis tem¬ 
plate, and we then total the x^ over all pulses. The results, 
shown in Fig.[^ indicate a dramatic improvement in x^ with 
only a few basis functions. In fact, the total x^ scales roughly 


as N, 


- 1/4 


though this relation is certainly coincidental and 


must steepen before Nt becomes equal to the total number 
of single pulses! 

It is implicit in our definition of a single-pulse template 
basis that each pulse corresponds to one and only one ba¬ 
sis function. We further make the assumption that pulses 
are uncorrelated. This is not true on short time-scales (a 
few pulses, Krishnamohan & Downs fl983 1, but is a good 
assumption for longer intervals. (We discuss pulse-to-pulse 
correlation and its mitigation further in !|^) With these as¬ 
sumptions, it is straightforward to simulate pulse trains and 
compare these with observations. First, we classify the 10"* 
single pulses according to basis function by minimizing x^- 
The fractions of pulses falling in each class are just the pa¬ 
rameters of a multinomial distribution. Second, we have ob¬ 
served in our sample that that the distribution of pulse flux 
densities within each class is approximately described by a 
lognormal distribution. Thus, we realize N simulated pulses 
by drawing an A-pulse sample from a multinomial distribu¬ 
tion and then, for each classified pulse, drawing a flux from 
the appropriate lognormal distribution. An example train of 
200 single pulses is shown alongside those taken from real 
data in Fig. 

Qualitatively, the simulations do a reasonable job of re¬ 
producing the observed pulse trains, though the real pulses 
show more stucture, particularly along the trailing edge of 
the peak. To quantify this, we have computed the modula¬ 
tion index (see |Jenet, Anderson fc Prince||2001[ for discus¬ 
sion), defined as m = [S) ~ /uiS)i i-®- the excess 


Figure 4. An example of 200 single pulses simulated as described 
in §4.2| Simulated pulses appear on the left, while 200 consecutive 
single pulses from PSR J0835—4510 appear to the right. 



Figure 5. The phase-resolved modulation index (see main text) 
of the flux density of 10^ single pulses from simulations (dashed) 
and data (solid). 


variance scaled to the mean intensity as a function of pulse 
phase. Figure shows the results for single pulses: most of 
the modulation in the leading peak is captured, while vari¬ 
ations are in the trailing region of the pulse in the real data 
are roughly twice those of the simulations. Much more in¬ 
terestingly, Figure shows a similar factor of two deficit in 
simulation variability for 100-pulse sub-integrations. While 
some of this discrepancy results from failing to include self¬ 
noise in the simulations, the dominant effect must be pulse- 
to-pulse correlation, e.g. as observed by [Krishnamohan ^ 


Downs (19831. On the other hand, we expect the simula¬ 


tions to become more faithful (approaching the single-pulse 
results) for longer sub-integrations. 


5 RESULTS 

Here, we describe the application of three different timing 
algorithms and the ‘pulse shape’ basis. The first algorithm, 
‘profiling’, does not require single pulse data, and we use it 
to introduce our analysis framework. 
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Figure 6. The phase-resolved modulation index of the flux den¬ 
sity of 100 100-pulse sub-integrations from simulations (dashed) 
and data (solid). 

5.1 Applications without single pulse data — ‘P’ 

In the most common mode of pulsar observation, the pul¬ 
sar signal time series is folded at the topocentric rotational 
frequency, producing a sub-integration containing many co¬ 
added single pulses. A typical integration time of 30 s records 
of order 100 pulses from a garden variety unrecycled pulsar. 
While the information contained in the single pulses is lost, 
we can still make inferences about the single pulse content of 
the sub-integration by estimating the coefficients of the ba¬ 
sis template. Those inferences in turn can furnish improved 
constraints on the phase offset S. In terms of the basis func¬ 
tion coefficients, the pdf for the sub-integration profile p 
becomes 

Pi-js- b)j.(^) y 
CTi ) 

(5) 

For pulsar timing, the basis coefficients s are nuisance pa¬ 
rameters. To maintain similarity with current pulsar timing 
techniques, we choose to use the profile likelihood to estimate 
5. That is, for each trial value of 5, we determine the value 
of s that maximizes the likelihood, reducing the likelihood 
to a one-dimensional function Cp{5). Because s is formally 
positive semidefinite (no negative flux densities!), we restrict 
the optimization of s to positive values. 

In the presentation below, we ignore the pulsar self¬ 
noise contribution and assume a constant a. We repeated 
our analysis including self-noise in the model, but found neg¬ 
ligible changes to the measurements of 5 and its uncertainty. 
That is, for Vela, the self-noise effect is small compared to 
pulse-to-pulse variability. As introducing the self-noise (het- 
eroscedasticity) substantially complicates the computation 
for little gain, we neglect it henceforth. 


5.1.1 Simulations 

To validate and characterize the approach under controlled 
but realistic conditions, we simulated pulse trains from 
the basis as described in §4.2| and then synthesized sub¬ 
integrations with 100, 200, and 400 pulses. These lengths. 



Figure 7. The precision of the phase measurement <5 as a func¬ 
tion of the complexity of the basis for three different realizations of 
the likelihood, using 100-pulse sub-integrations of simulated data. 
The grey, solid trace (‘P’) corresponds to the profile method, dis¬ 
cussed in The dashed salmon and dot-dashed black lines cor¬ 
respond to the ‘tagging’ methods discussed in §5.2| both without 
(‘T’) and with (‘M’) the addition of the multinomial likelihood. 
All results are scaled to the ideal measurement precision (see main 
text). The inset is as the main figure, but zoomed to show that 
the tagging methods perform optimally, reaching the Fisher pre¬ 
cision. In these results, the abscissa indicates for both the 
simulations and the fitting, i.e. as if one had perfect knowledge 
of the family of single pulse variations. Finally, note that the er¬ 
ror bars due to finite sample size are shown for the Nb = 1 and 
suppressed for the remainder; see main text for details. 



Figure 8. As Fig. but the simulations use an Nb = 100 basis 
while the abscissa indicates the Nb used in the fits. This case cor¬ 
responds to imperfect knowledge of the panoply of pulse shapes. 


and the number of simulated sub-integrations, are chosen to 
match those of the data, described below. 

For each sub-integration, we maximized the likelihood 
of equation]^ to determine the phase shift estimator d. We 
then computed the standard deviation of S over the set of 
simulations and normalized it to the idealized scatter, i.e. 
what we would observe in the absence of jitter/SWIMS and 
self-noise. This uncertainty is given by the Fisher informa¬ 
tion of the mean template, which gives a lower limit for the 


Nb 


/(p I (5, S, b, cr) = P 


7 = 1 ^ 


■ exp 


■7T(Ji 
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variance for 5, 


-2 _ / 1 dti 

'^,5,Fisher ^ I Qg 

i ^ 

where Cr is again the radiometer noise. 

In our simulations, we can choose to simulate profiles 
with a basis with the same Ni, we use to maximize the likeli¬ 
hood, or we can fix it at some representative Nb- The former 
essentially validates the method, while the latter probes how 
the results change when the basis and the data disagree, as 
will certainly be the case for real data. Accordingly, we fol¬ 
low both approaches. 

The results for which Nb is the same for both simula¬ 
tions and fitting appear in the grey traces labelled ‘P’ in Fig. 

Here, for clarity, we show the results for the 100-pulse sub¬ 
integrations; the 200-pulse and 400-pulse results are similar. 
The leftmost values, corresponding to a template basis with 
Nb = 1, i.e. the mean template or ‘traditional’ pulsar timing, 
shows a scatter consistent with the Fisher limit—certainly to 
be expected, since simulations from a template with Nb = 1 
have no jitter/SWIMS. Since the set of simulations changes 
for each Nb, the data are uncorrelated; this will not be the 
case for scenarios described below. As the complexity of the 
simulations increases, we see the profile method slowly di¬ 
verges as the degrees of freedom increase. Thus, we expect 
the method may ultimately break down if we attempt to fit 
a pulsar with very complicated pulse-to-pulse variability. 

In the second set of simulations, we fix Nb = 100 for a 
high level of pulse-to-pulse variability, and we then vary Nb 
in our fits of 5. Now, the Nb = 1 case shows a jitter/SWIMS 
level of about 50 times the radiometer limit, consistent with 
our observation that our simulations produce about half of 
the pulse-to-pulse variability of the actual data. As we in¬ 
crease Nb, the profile method is able to successfully recon¬ 
struct some of the single pulse variability, reducing the jit¬ 
ter/SWIMS level by a factor of three for Nb = 20, before 
slowly diverging. Since we use the same set of simulations for 
each value of Nb on the abscissa, these points are correlated, 
and the error bars should be interpreted as uncertainty on 
overall amplitude rather than as a measure of scatter. 

Thus, we conclude that this ‘profile’ method is quite 
effective in removing jitter/SWIMS noise under the ideal¬ 
ized conditions of our assumptions about single pulse vari¬ 
ability. This approach should be especially effective for pul¬ 
sars with less complicated pulse-to-pulse variability, and for 
those whose single pulses are too faint for effective single 
pulse ‘tagging’ (see below), such as MSPs. 



5.1.2 Data 

We now apply the profile method to synthesized sub¬ 
integrations of varying lengths from the data set described 
above. After excising the pulses used to build the basis 
template and those affected by RFI, 11,182 pulses remain. 
Maintaining a sufficiently large sample of sub-integrations to 
measure the scatter of 5 imposes a limit on sub-integration 
length of about 400 pulses, or 36 s, typical of a fold-mode 
sub-integration. However, as we discuss below, we expect 
the results to scale to longer integrations, and perhaps to 
improve as correlations between single pulses become less 
important. 



Figure 9. As Fig. but for data from PSR J0835—4510. Addi¬ 
tionally, the orchid, dotted trace shows the results for the flux- 
tagging scheme described in §5.3| 


As with the simulations, for each synthesized sub¬ 
integration, we maximized £.p{S) to determine 5 and exam¬ 
ined the distribution of values. The results are summarized 
in Fig. As before, the leftmost points correspond to ‘tra¬ 
ditional’ pulsar timing and show the jitter/SWIMS level is 
about 100 times the ideal limit. As with the second set of 
simulations, because the same data set is used for each set 
of basis templates, the results are correlated. 

As we might expect from the failure of our basis to re¬ 
produce the lion’s share of the pulse-to-pulse variability, the 
reduction in scatter in moving to Aij, > 1 is substantially 
poorer than in simulations. Nonetheless, we see a gradual 
impovement in the level of scatter, as the basis complexity 
increases and captures more of the pulse-to-pulse variabil¬ 
ity, with a typical scatter reduction of about 25%. We also 
note that longer sub-integrations show a lower level of jit¬ 
ter/SWIMS relative to the ideal (which scales as l/\/t) than 
shorter sub-integrations, suggesting that performance does 
indeed improve as the pulse-to-pulse correlation time-scale 
becomes short compared to the sub-integration length. 

While a reduction in the scatter of measured phase off¬ 
sets clearly improves precision, an accurate estimate of the 
remaining noise is of vital importance for pulsar timing. Ide¬ 
ally, reliable error estimates can be extracted directly from 
the timing procedure, and we expect this could be achieved 
even without a perfect template basis via bootstrap meth¬ 
ods. But minimally, all we require is that the normalized 
residuals follow a gaussian distribution. These normalized 
residuals appear in Fig. |10| showing that the estimators for 
S do indeed follow a normal distribution. 

5.2 Applications with single pulse data — ‘T’ and 
‘M’ 


Next, we consider a hybrid mode of observation—one in 
which we produce sub-integrations as before, but in which 
we also have access to single pulse data. While many avenues 
are available, including marginalization over template classi¬ 
fication and Monte Carlo Markov chain methods, we opt for 
a simple approach that might more readily be implemented 
in an online pulsar data recorder. Because nearly every pulse 
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Figure 10. The distribution of the error-normalized residuals af¬ 
ter normalization by the sample standard deviation, for all three 
‘shape’ algorithms, with 100-pulse sub-integrations and Nf, = 20. 
The black trace and error bars show the expected gaussian dis¬ 
tribution and the range covered by 90% of random realizations 
of identical histograms of gaussian random variables. The colors 
follow the convention of Fig. The orchid line shows the same 
result for the ‘flux basis’ of |5.3| 


in our data set has a high S/N, it can be uniquely associ¬ 
ated with a single basis function. This leads to the idea of 
‘pulse tagging’—each incoming single pulse is classified and 
its flux (scale parameter) measured; the pulse itself can then 
be discarded. The classihcations and scales are then used to 
build an optimized template from the template basis, which 
is in turn used to measure the phase offset of the full sub¬ 
integration. A key point: we are not “using the data to ht 
the data”. Because we do not know the correct phase offset 
S a priori, we perform classification and template construc¬ 
tion over a range of values of S. This process is similar to 
constructing a profile likelihood, in which we optimize for 
the (discrete) profile classification parameters. 


5. 2.1 Simulations 

As before, we test this method with simulated data. In the 
hrst ‘sanity check’ approach, we see that the approach per¬ 
forms flawlessly. That is, by tagging the individual pulses, 
we are able to construct a template for each sub-integration 
that entirely accounts for jitter/SWIMS, and we reach the 
Fisher limit for estimator uncertainty. Importantly, despite 
implicitly fitting for the proHle class and scale of each pulse, 
we see no evidence for overhtting. 

The second set of simulations, in which all simulated 
profiles are drawn from a Nf, = 100 template basis, are more 
interesting. For hts with Nf, 100, there is no improve¬ 
ment in jitter/SWIMS level. But once we reach Nf, = 50, 
jitter/SWIMS is reduced by roughly half, and as in the hrst 
set of simulations we remove jitter/SWIMS entirely once Nb 
is consistent between simulation and ht. This result is read¬ 
ily explained. For the Nb 100 cases, no template basis 
function maps well onto a single pulse, and tagging process 
will essentially build up the mean template. On the other 
hand, once the basis begins to encapsulate single pulses, we 


achieve a template that accurately records the inhuence of 
pulse-to-pulse variability. 

From these results, we conclude that the single pulse 
tagging approach is extremely ehective (entirely removing 
jitter/SWIMS) if the template basis is well-suited to the 
data. 

5.2.2 Data 

The results of applying this approach to our test data set ap¬ 
pear in the salmon traces of Fig.j^ Jitter/SWIMS decreases 
rapidly with increasing Nb through Nb = 5, before saturat¬ 
ing at roughly 65% of the Nb = 1 level. These results are 
largely consistent with our understanding that the template 
basis should be well matched to single pulse variability. As 
we saw in Fig.[^ we capture most of the pulse-to-pulse vari¬ 
ability that we can with relatively few template basis func¬ 
tions, and increasing Nb beyond about 5 yields only slow 
gain. 

Although the tagging method seems effective in simu¬ 
lations, there is in principle additional information we can 
apply to the fits that may better handle real data. E.g., very 
rare, bright pulses can substantially affect a TOA measure¬ 
ment ( [Oslowski et al.|2014p , so they must be represented in 
an effective template basis. On the other, less intense pulses 
with morphology similar to giant pulses should not be so 
classified. Put another way, the pulse counts and flux den¬ 
sities associated with each basis are random variables, and 
we can add terms reflecting this to the likelihood. Substan¬ 
tial misclassification will then be penalized. As discussed in 
§4.1| we can easily determine the expected rates of pulses 
associated with each basis, as well as measure their flux 
distribution. With our neglect of pulse-to-pulse correlation, 
these terms enter the likelihood multiplicatively through the 
multinomial distribution (for the classes) and Np lognormal 
distributions (for the fluxes). Suppressing constant terms, 
this log likelihood is 

1 fl 

\ogC = --'^[pi-U{S)f + - 

i=l 

Nb 

+ - log rii!], 

i = l 

where 

Np 

t = ^Sibi, (8) 

i = l 

bi is the best-ht basis for the ith pulse, and Si is the cor¬ 
responding best-fit flux. The a; are the parameters of the 
multinomial distribution (relative pulse frequency for each 
basis) and the Ui are the number of pulses tagged into each 
basis. 

In our data—and in any data in which jitter/SWIMS 
is important—the ‘x^’ term in equation utterly domi¬ 
nates the multinomial and lognormal terms, and the relative 
weighting must be adjusted to allow the latter data to pro¬ 
vide constraints. There are two ‘natural’ approaches, and we 
Hnd they yield similar results: (1) re-scale the term such 
that the residual per degree of freedom is unity or (2) 
re-scale the x^ such that the x^ term contributes equally to 
the total likelihood. We adopt the latter approach. 


Dg Si — fib 


CTb 


+ 2 log Si 


(7) 
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Applying these additional constraints yields a small de¬ 
crease in the scatter of our sample, in some cases reducing 
the jitter noise by as much as 40%; see the black traces of 

Fig-H 


5.3 A Flux Basis 


The previous methods have all focussed on pulse shape 
rather than information encoded in the pulse profile normal¬ 
ization, i.e. its flux density. However, for PSR J0835—4510, 
there is a additionally a strong correlation between the peak 
intensity of a pulse and the phase at which that peak falls 


(Krishnamohan & Downs 19831, and strong pulses have a 


large influence on TOA estimation. We thus consider here 
a ‘flux basis’ and investigate its performance relative to the 
‘shape basis’. 

We construct this basis by simply computing the peak 
S/N in each single pulse, dividing the pulses into bins loga¬ 
rithmically uniform in peak S/N, and constructing templates 
from the average pulse shape in each bin. To perform tim¬ 
ing on a synthesized sub-integration, we assume as in §5.2| 
that we can measure single pulse properties, specifically the 
peak S/N in each single pulse, selecting the basis template 
from the matching S/N bin, and adding up the matched tem¬ 
plates. The resulting template naturally reflects the outsized 
contribution of strong fluxes to a sub-integration. Note that 
it also establishes the additional observational requirement 
of stable gain/flux calibration. 

The results are extremely interesting. In Figure we 
note that the flux basis captures less of the pulse-to-pulse 
variability. The timing performance, however, is excellent. 
The results appear as the orchid-coloured trace labelled ‘F’ 
in Figure showing a performance that is generally the 
best of the methods explored here, offering a robust 40% 
improvement over the jitter/SWIMS noise level! 

The root of this discrepancy between superiour timing 
performance and poor pulse-shape reproduction offers some 
insight into the Vela’s pulse emission process. Aside from the 
clear flux/phase correlation, it is evident that the emission 
mechanism can produce pulses with similar shapes at differ¬ 
ent observed phases, which in turn may point to propagation 
effects as an important source of jitter/SWIMS. Whether 
similar effects are common in other pulsars, and how they 
might be explored (e.g. through polarimetry) are interesting 
topics for future work and we discuss them further below. 


6 SUMMARY AND FUTURE WORK 

We have presented a new method for pulsar timing us¬ 
ing single-pulse based likelihood. In simulations, we showed 
these methods are capable of materially reducing or en¬ 
tirely eliminating the jitter/SWIMS noise associated with 
pulse-to-pulse variability. We demonstrated less dramatic 
but nonetheless substantial reductions in the jitter/SWIM 
noise of data from the bright Vela pulsar. 

Our results are encouraging and motivate future work. 
In particular, much could be done simply be expanding the 
scope of the data used to build the basis function. Our ‘train¬ 
ing set’ of pulses for Vela was limited to a mere single 
pulses, comparable in size to the set of pulses we could use 
to test the algorithms. It must necessarily have missed rare 


pulses, and the limited number of basis function likely forces 
the combination of similar but distinct pulse classes into a 
single template. Indeed, our simulations, in which we had 
perfect control over the complexity of pulse-to-pulse vari¬ 
ability, showed that we required a basis with nearly the full 
complexity of the simulations to fully remove jitter/SWIMS 
(Figure [^. 

Moreover, Vela exhibits pulse-to-pulse correlation, con¬ 
trary our simplifying assumption of independent pulses. 
Though more complicated computationally, a ‘two-pulse’ 
template basis—whose members are two consecutive single 
pulses—is straightforward conceptually and could capture 
much of the observed single pulse correlation. Extension to 
‘N-pulse’ bases would help even more. Both cases would re¬ 
quire more extensive training sets. Likewise, it is straight¬ 
forward to extend these methods to pulsars with periodic 
sub-pulse modulation, in which case the natural unit of a 
basis is the modulation period. 

In general, it seems likely that per-pulsar tuning will 
be required to find the best combination of template ba¬ 
sis and fitting method. For Vela, an optimal method must 
clearly take advantage of the TOA / flux density correla¬ 
tion. In the case of PSR J0437—4715, Oslowski et al. (20131 
observed a correlation between the phase centroid and the 
pulse polarization, suggesting an optimal template basis for 
thise pulsar must include polarimetry. For mode changing 
pulsars, a separate template basis for each mode could be 
determined, and subsequent application of the timing algo¬ 
rithms would furnish both correct TOAs and an automatic 
classification of the pulsar modes. 

While there is much future work to be done, a press¬ 
ing question remains: if the current template basis can do a 
good job of describing single-pulse variation (Figure]^, why 
is the timing only modestly improved? This ultimately boils 
down to the still poorly understood phenomenology of single 
pulse variability. For pulsars whose emission truly comprises 
‘shots’ of emission within a guidin g envlope ||Cordes||1975p , 
fully modelling the covariance a la Oslowski et al. (20111 is 
the optimal approach. Whereas pulsars which demonstrate 
at least some underlying pattern, in which similar single 
pulses appear repeatedly over time, will benefit from the 
recovery of additional information by the deterministic ap¬ 
proach advocated here. Though the pulse-to-pulse correla¬ 
tion observed for PSR J0835—4510 implies at least some 
repetition, it may yet be that the single pulses also possess 
an irreducible random component. In this case, the reduc¬ 
tion of jitter/SWIMS by 30-40% achieved here may be the 
best possible, and the rest that can be done is accurately 
estimate uncertainties. Indeed, improving the timing algo¬ 
rithms in the ways outlined above will provide a stringent 
test of the true randomness of the pulsar mechanism, and it 
is thus of keen interest to collect single pulses from a large 
set of bright pulsars for further study. 

What would be required to implement this method in a 
practical pulsar timing programme? As we have discussed, 
a sina qua non is a large training set of high-quality single 
pulses. These can be obtained in a single dedicated session, 
possibly on a target-of-opportunity basis when a pulsar is 
unusually bright due to scintillation. Subsequent observa¬ 
tions should record single pulses, but possibly with reduced 
quality. Compared to typical folded sub-integrations, full- 
resolution single pulse data from young (millisecond) pul- 
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sars requires about 100 (10,000) times more storage. While 
full-resolution data is intrinsically interesting for studying 
the pulsar mechanism, our method is based on fully re¬ 
duced (in frequency and polarization) data, and recording 
such archives requires only modest additional storage. We 
strongly advocate such a commensal programme, as it will 
both enable our proposed study of jitter/SWIMS and help 
inform the design of future pulsar timing instrumentation. 
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